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ABSTRACT 


An amphora catalyst which resembles a hollow sphere 
with two openings, on the two sides, has been cho sen for study. 

The active catalyst is dispersed nonuniformly over amphora. Two 
types of distribution have been studied, the first one is non - 
uniform distribution with a, as activity distribution parameter 
whose value has been changed from 0 to 8 and the second one Dirac- 
delta type distribution. 

It has been found that amphora catalyst with uniform 
activity distribution (a=0) gives better effectiveness factor and 
selectivity over conventional spherical catalyst. When the value 
of a is increased to 1 (linear distribution), or to 2 (parabolic 
distribution) or more than 2, the activity and selectivity is 
found to increase. 

*For Dirac-delta distribution, generalised optimal 
location has been defined and it is shown that as the thiele 
parameter increases, the optimal location shifts towards the 
surface of the pellet and at very high values of thiele parameter, 
the optimal strategy is to apply all the active material on the 
surface. 
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NOMENCLATURE 

Boundary value at the enner surface. 

Volume averaged catalyst density. 

Local catalyst density. 

Boundary value at the outer surface. 

Concentration of componentsA and B respectively. 
Surface concentration of componentsA and B, respective 
Effective diffusivities of componentsA'iand B rpspecti\( 
A vector containing the boundary values at external 
surface* 

A function of u^. 

Step width. j 

Local rate constant for first and second reactions. 
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Global selectivity. | 
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Initial surface concentration of component B. 
Dimensionless concentration gradient of components 
A and B respectively. 

Concentration of component A and B respectively where 
the active catalyst is applied. 

Optimal concentration of component A. 

Volume element. 

Volume of the pellet. 


V 



(iv) 


Y 

z 

^opt 

a 

P 

d 

^ext»^int 


01,02 

Q 

a 

A 


Vector containing boundary values. 
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INTRODUCTION 


Mass and heat transfer limitations are often 
important considerations in heterogeneous catalysis. Chemical 
reactions that are diffusion limited are being greatly benifi- 
tted from the use of smaller catalyst particles, because of 
their smaller internal diffusional path length. By minimizing 
residance time of product molecules within the catalyst parti- 
cles, secondary reactions can ‘be minimized and selectivity can 
be increased. 

Keeping in view, the higher surface to volume ratio 
the researchers have introduced a variety of catalyst shapes 
over the years. Spherical catalysts have been used in place 
of extrudates or tablets because of their higher surface to 
volume ratio and better packing efficiency. 

THE AMPHORA 

In the late 1950* s Sohio research laboratory was 
busily expanding its acrylonitrile technology, in course of 
developing a spray drying process for making fluid bed 
catalyst (17), 

Jim Callahan made an interesting discovery, upon 
viewing a sample of spray dried catalyst under the microscope. 
He found that a number of particles had an unexpected shape. 
Scattered among the spheres were particles that had channels 
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Amphora I 



fig. lA Cross - sectional drawing of a representative 
particle amphora I 


rig. IB Cross-sectional drawing 6l a representative 
particle amphora 11 



Fig. 2 -Mechanism of amphora formation. 





leading from the outside to the hollow centre. These new 
particles, if they can be produced in sizes appropriate for 
fixed-bed applications, would offer a unique combination of 
advantages. The new shapes have a higher surface to volume 
ratio, than solid spheres . Most importantly, if the inner 
cavity would be made to be half the diameter of the particle, 
the effective particle thickness would be l/4th that of the 
particle diameter. Thus a bed of 3 mm diameter particles 
would have the performance characteristics of 0.75 mm solid 
sphere, Fig. (1). 

Name for the new particle shape, again, was sugg- 
ested by Callahan, He called them Amphorae, for the ancient 
Greek urns that were used to carry cargo on trading ships. The 
Greek amphorae were shaped so that they stacked uniformly and 
efficiently in the ships holds, and therein lies the analogy. 

The industrial process of making amphora was described by 
Harrington (17). 

MECHANISM OF AMPHORA FORMATION 

The actual mechanism for amphora formation, develo- 
ped by G.K. Meloy, is a variation of one of the classic mechanism 
of the drying of a sphere in spray drying processes. In that 
mechanism a rigid shell is initially formed on the surface of the 
sphere so that further drying occurs by migration, solid parti- 
cles are carried with the liquid. Because the rate of liquid 
migration and evaporation is faster than the rate of diffusion 
of the solid particles back into the interior, a hollow sphere 
is formed. * 
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Meloy found that the variation of this mechanism 
that results in amphora formation is one of unequal evapora- 
tion rates over the surface of the spherical particle. This 
is attained by a unidirectional air flow over slurry droplets 
as they are drying. The air flow speeds evaporation on one 
side of the exposed portion of the droplet, rigidizing that 
portion but leaving an area of greater fluidity on the other- 
side of the droplet capable of plastic flow. As drying conti- 
nues and natural shrinkage occurs, the particle first forms a 
hollow centre and then, to relieve surface stresses and mini- 
mize surface energy, the area of greater fluidity shrinks into 
the particle and the hole is pulled. Fig. (2) . 

Experiments carried out at Sohio (17) proved that 
amphora catalysts could indeed outperform conventionally shaped 
catalysts in reactions limited by heat and / or mass transfer. 
The tables 1, 2, 3, compare the amphora and other conventional 
catalysts. 

INTRODUCTION TO NON UNIFORMLY ACTIVE CATALYST 

An extensive literature exists dealing with the 
interactions between the physical and chemical transport 
phenomena of the heterogeneous catalysis and the distribution 
of catalytically active material within the particles assumed 
to be uniform ( 3,4,31,44,45). The effects of intraparticle 
mass and heat diffusional resistances in catalysts exhibiting 
uniform activity distributions have been extensively studied by 
these workers in terms of point and overall effectiveness. 
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-Various typ^s of distributions. 
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selectivity, yield and the concentrations and temperature 
profiles, within a single particle (34) and in the chemical 
reactors. 

Kasaoka and Sakata (22) first studied the non 
uniform catalytic effectiveness factor with variable rate 
constants and diffusiv ities theori tically on flat pellets. 

Corbett and Luss (8) have shown the effect of non uniform 
activity distribution of a single spherical pellet on its 
effectiveness factor, selectivity and deactivation • In most 
cases the performance of nonuniformly active catalysts is 
superior to that of uniformly active catalysts which have the 
same volume averaged activity. The specific deactivation 
mechanisms and rate as well as diffusional resistance have an 
important influence on the most desired activity profile. 

Bimolecular Langmuir kinetics under isothermal 
conditions, was studied by Becker and Wei (5) and shown that 
the *egg yolk' type distribution gives the greater activity 
than that the traditional egg shell or uniform distribution 
Fig. (3). The optimal activity distribution inside the 
isothermal and nonisothermal pellets of spherical, cylindrical 
and slab geometries with negligible and considerable external 
mass transfer resistance under steady state conditions was 
studied by Morbidelli et.al (24,25), Juang (19) Ready Karri (30), 
Gottifredi (14), Shadman-Yazdi (32), and Vayenas and Pavlov (35,36, 
37). These researchers concluded that the optimum distribution 
is the approximately chosen Dirac -d type distribution and they 
generalized the optimum distance for varying conditions, giving 
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the maximum effectiveness factor, selectivity, yield and 
durability of the catalysts. Analytical expressions for 
this location as the function of all the involved physio- 
chemical parameters are given in three pellet geometries 
by Vayenas and Pavlov (35,36,37)* Numerical computations 
are reported by Morbidelli . Servida (24) to analyze the case 
where the active catalyst is deposited in a step manner center- 
ed around the optimal location. 

Same results were ' obtained by Morbidelli (27) with 
reference to the synthesis of ethylene oxide. The most suitable 
objective function has been proposed and a very easy optimization 
technique, which identify a close approximation of the optimal 
radial activity distribution, has been developed. This leads to 
a Dirac del ta , function in practice approximated by a step distri- 
bution whose width is about 4yi of the particle radius, centered 
at a particular radial location within the particle, which depe- 
nds upon the operating conditions. 

Johnson and Verykios (20) studied theoritically the 
selectivity enhancement in ethylene oxidation employing partially 
imprenated catalyst, using the experimental data and reported the 
effects of radially nonuniform distributions of catalytic activity 
on the performance of spherical catalyst pellets under ethylene 
oxidation. They also showed that partially impregneated cataly- 
st exhibits higher activity and selectivity at higher value of 
thiele modulus, than the uniformly active catalyst pellets. 'They 
also observed the multiple steady state effectiveness factor 
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and selectivity under sufficiently high values of the heat 
of reaction parameter. 

Simultaneous, bimolecular, complex reactions, 
and Michaelis -Menton kinetics under exothermic conditions 
was studied by Johnson and Verykios (21) and Cukierman et.al (9), 
Horvath (18) . They assumed the distribution function of the 
form: k(r) =ks(r/r-)°^, a>0 and the reactions were assumed to 

o 

be of the second order. Effectiveness factor and overall 
pellet selectivities for various values of the dimensionless 
activation energies and heat of reaction parameter were repor- 
ted as a function of the thiele modulus. 

Analytical evaluation of the isothermal and non- 
isothermal fixed bed reactors were reported by Morbidelli and 
Verma (28,29), Smith (33),Juang and Weng (19), Cukierman et.al (9). 
They showed that the reactor can possess at most three steady 
states and that the issue of reactor multiplicity and stability 
is intimately tied up with multiplicity and stability for the 
first layer of pellets packed along the axial direction. The 
patterns of reactor behaviour, including ignition phenomena is 
fully characterized as a function of three dimensionless parame- 
ters viz, adsorption constant, p, which depends upon thiele 
parameter and the Damkohlar number. 

For a bimolecular Langmuir-Hinshelwood kinetics, the 
rigorous optimal distribution is found to be a Dirac delta func- 
tion centered at a location inside the pellet which depends 
on position along the reactor axis. They have also shown that 
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the nonunif ormly active distribution also plays an important 
role in controlling the temperature profiles inside-the reactor. 
Smith (33) reported that yield can be substantially improved 
by a step distribution. 

Poisoning with shell progressive models, poremouth 
poisoning, activation policy for poisoning, parallel and series 
poisoning was studied by Vayenas (38), Ardiles (2), DeLancey (10) 
and Juang and Weng (19). Vayenas (38) reported that when the 
rate limiting step is the chemical reaction the optimal strategy 
is to concentrate the active component in a thin shell located 
inside a pellet. Thus the poisoning by shell progression can 
be avoided. The same results were reported by Ardiles. Juang 
and Weng shown that the convex and concave profiles may signi- 
ficantly resist the homogenous poisoning in the form of series 
and parallel fouling. 

Shape normalization and the analytical expressions 
for the effectiveness factor and selectivity and their analytical 
solution were given by Yortsos and Tsotsis (42,43), Wang and 
Varma (39,40), Morbidelli and Varma (26), They developed the 
shape normalization technique for nonuniformly active catalyst 
pellets with an isothermal first order reaction. The shape 
normalized curves, obtained for different geometries, activity 
profiles and external mass transfer resistance match perfectly 
for small and/or large values of the thiele modulus. Yortsos 
and Tsotsis developed the shape normalization by using the 
Li ouvill e-Green expression for asymptotic behaviour of the 
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effectiveness factor. They also reported that the normaliz- 
ation has no meaning when the surface catalytic activity is 
zero. 

The optimization by very simple technique was done 
by Morbidelli (24,25), DeLancy (10). They employed the 
Pontryagan's continuous maximum principle to obtain the optimal 
profiles. Dougherty and Verykios (12) reproduced the various 
techniques of optimization and reported that a^Dprximately chosen 
Dirac delta is the best recommended profile. 

Beyond the theoritical research impregnation models 
and experimental studies, on the nonuniformly active catalyst 
was done by Harriott (16), Minhas and Carberry (23), and Wu 
et.al (41). Wu, first time studied experimentally, the exis- 
tence of the optimal active catalyst distribution, and the 
dependence of the effectiveness factor vs. thiele modulus in 
nonisothermal pellets. For the hydrogenation of ethylene on 
Pd/Al202 catalyst and the methanation of Co on Ni/Al202 catalyst, 
by use of step ditribution.with the active layer of narrow width 
to approximately replace the distribution of the Driac delta 
type, it was shown that there did exist an optimal location of 
the active layer within the pellet at the given conditions and 
the optimal location moves towards the external surface as the 
thiele modulus increases. 

In the present work we theoritically studied the 
active profiles in the new catalyst shape, called amphora. 
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We reported the effectiveness factor, selectivity and yield 
for the various type of distributions like uniform, linear, 
parabolic, hyperbolic and the Dirac delta. We have shown that 
the maximum selectivity and the effectiveness factor was obta- 
ined by the Dirac delta distribution, the site of which depends 
upon the operating conditions. It has, also been shown that 
at very low value of thiele modulus the whole thing should be 
centred in the pellet, and as the thiele modulus increases 
the optimal location moves towards the surface and thus when 
the reaction is highly diffusion controlled, i.e., when the 
thiele modulus is very high, the whole active sites should 
be applied on the surface. Which is properly known as the 'egg 
shell' type distribution. 

A generalized activity distribution function, has 
also been defined. byhequati on (25) . Here the most important 
thing is that, since the amphora has the two surfaces, hence 
we need two type of profiles in the same pellet. And also in case 
of the Dirac delta type distribution there exist two' optimal 
locations, which has been reported. 
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Table 1. Oxidation or propylene to acrolein/acryllc acid 
(amphorae vs tableted catalyst) 




Corrected ^ per 

pass conversion 

Catalyst form 

Temperature 

(^C) 

Acrylic 

Acid 

Acrolein 

Total ; 

Conver- 
sion 

Useful 

Smin X 2.8 mm 
Tablets 

310 

4.8 

74.8 

81.7 

79.6 

5mm X 2.8 mm 
Tablets 

329 

7.7 

82.6 

92.9 

90.3 

5mm X 2.8 mm 
Tablets 

338 

8.8 

82.0 

93.6 

90.8 

3.1mm Diameter 
amphora 

310 

7.1 

84.3 

93.9 

91.4 

3.1mm Diameter 
amphora 

327 

10.8 

83.1 

96.8 

93.9 

3.1mm Diameter 
amphora 

338 

13.1 

82.6 

99.5 

95.7 
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Table 2. Oxidation of Oxylene to phthallc anhydride 
(amphorae vs. solid spheres) 





Corrected '/ per pass 

Conversion 

Catalyst 

Temperature 

Exotherm 

Phthallic 

CO 

co^ 

from 

(^^c) 

(°c) 

anhydride 



Amphora 

353 

18 

67.8 

3.3 

11.4 

Spherical 

348 

20 

71.2 

7.3 

21.5 

Amphora 

368 

•20 

78.8 

5.3 

14.1 

Spherical 

361 

26 

63.8 

' 9. 5 

26.7 

Amphora 

381 

22 

77.7 

5,9 

16.4 

Spherical 

376 

40 

53.3 

12.9 

33.9 

Table 3. 

Oxydehydrogenation 

of 2-butene to butadiene (amphorae 


vs. tablets) 





Catalyst 

Temperature 

Butadiene Maleic 

Total 

useful 

form 

(°c) 


anydride 

conversion 

Tablet 

302 

40.4 

1.1 

44.5 

41.5 ; 

Amphora 

302 

48.2 

1.5 

53.5 

49.7 J 

Tablet 

329 

69. 5 

3.4 

78.5 

72.9 

Amphora 

329 

72.6 

6.7 

88.0 

79.3 


Reported from reference No, (17), 
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FORMULATION OF THE PROBLEM 


As it is seen the Amphora II catalyst Fig.(l) 
particles are spherical in shape and the most advisible ratio 
of inner to outer radius is 1/2, as was suggested by Herrington (17). 
Amphora II has two naturally occuring openings on diametrically 
opposite points. Keeping these features in mind about the shape 
of the Amphora II, following assumptions are made in modelling 
of the non uniformly active catalyst for diffusion problems with 
negligible external phase mass transfer resistance; 

1. The catalyst particle is spherical in shape with cavity 
inside. 

2. There are two openings for fluid transportation to and from 
the inside surface. 

3. Ratio of the external diameter to the internal diameter is 2. 

4. The volume occupied by the openings is small as compared with 
the total volume of the particle. Thus the small material loss 
can be ignored. 

5. The external phase mass transfer resistance, out side the 
external surface and in the cavity is negligible. 

On the basis of the above assumptions the model reduces 
to the hollow sphere with ratio of outer to inner diameter being 2. 
And the boundary conditions, on the basis of the assumption 5, are 
reduced to the Direchlet type. 

CONSECUTIVE REACTIONS: 

We consider the following first order, isothermal, 
consecutive reaction with negligible external phase mass transfer 
resistance. 
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A 



B 


K, 


» C 


Where and K 2 are the rate constants and B, the 
intermediate product, is the desired product. also , The steady 
state mass balance equations for A and B follows from Aris (3), 


Dej^V^Cj^ = a(r) (1) 

= {K 2 C 2 -K^Cj^) a(r) (2) 


Where Cj^ and C 2 are the concentrations of components 
A and B respectively and a(r) is the local active catalyst density, 
which is defined as 


/ a(r)dv = a^v (3) 

vp 

Where a^ is the volume averaged active catalyst distri- 
bution inside the pellet and v^ is the total volume of the amphora 
pellet. 

On the basis of the assumptions made in the beginning of 
this chapter, the boundary conditions, to solve the steady state 
mass balance equations (1) and (2) are given as 


^1 ~ *^10 surface (4) 

^2 “ ^20 surface (5) 


NONDIMENSIONALIZATION OF THE EQUATIONS ; 

We introduce the following nondimensional quantities to 
solve the equations (1) and (2) numerically. 
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^2 " 


— 2 


, 




^ .2 


I- r2 

24 ^ 


K, a(r) , K, 


^2 


a( 2 ) 


C ’ 

^10 

a(r) 


As the amphora II under certain simplified assumptions 
becomes the hollow sphere. Putting these dimensionless variables 
into equations 1 to 5 and expending the equations (1) and (2) in 
the spherical geometry, we got the following dimensionless 
equations. 


1 d 


2 du^ 


dT = ’^l “1 


) = 07 u, aCz) 


_2 dz dz ' 


02 U 2 a(z) - *^0^ Uj^ a ( 2 ) ... (7) 


and the boundary conditions reduces to. 


u^ = 1.0 on the surface 


U2 = U2 


on the surface 
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As was suggested by Herrington (17) that th'e most 
probable ratio of external to internal diameter is 2 , so the 
equations ( 8 ) and (9) can be written in the more general form. 


u^(l) = u^{0.5) = 1.0 ( 10 ) 

U 2 (l) = 02(0.5) = U2 ( 11 ) 


Equation (3) can be reduced in the dimensionless from, 
and each type of profile should satisfy this equation. 


1.0 >-j 

S aiz) dz = ^ 

0.5 


( 12 ) 


From equation (7) the nondimensionless rate of formation 
of component B, is given by, 

dz = ^ 

0.5 


i _ 

I Z dz = / (02 u^«^02 ^2 .,.(13) 


0.5 


It can be further simplified as , 
du ^ 


(^) -0.25 = / (02 ^ ^^02 j ^2 ..(14) 

2=1 O.O 


du. 


dui 


But ( 3 ^ ^z=l ~ ^ ^0 5* both are on the opposite 

directions and ofcourse at the surfaces. 

Hence equation (14) can be reduced 


^'^2n 4 0m2 .. ^2 


= 5 / (0f ^ 02 U 2 )a(z)z'^ dz ... (15) 

0.5 • 



18 


Similarly the dimensionless rate of consumption of component 
A is given by the following equation, 


du, 

(dT ) 


^ / 0^ u, a(z) 2 ^ d 2 

^0.5 ^ 


( 16 ) 


Selectivity can be difined as the net rate of formation 
of the desired product to that of the rate of consumption of 
the key component. In the mathematical form, the global selec- 
tivity of this single pellet can be defined as, 


S = - 


De,(d^) 

on surface 



^®1 ^ dr ^ on surface 


(17) 


Making this equation dimensionless and putting the 
values of the derivatives from equations (15) and (16) we got 
the following relations. 




du^ 

on surface 

du^' 

^ dz ^ on surface 


(18) 


, I 0o ® 

S = 1-1 Q. 5 ^ ^ 

/ 0^ u, a(z) z^ dz 

0.5 ^ 


(19) 


( 20 ) 


We further reduce equation (19) as following 


/ u« a(2) dz 

S = 1.0 - 0.5 ^ 

1 / X 2 

J u^ a(z) z dz 

0.5 


( 21 ) 
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The ratio of surface fluxes expressed as the 
slectivity S will be negative when there is the net consumption 
of the desired product and positive when there is the net for- 
mation of this product* 


The effectiveness factor of the component A is defined 
as ratio of the rate of diffusion effected reaction to that of 
diffusion free reaction. More mathematically this can be 
represented as. 




1.0 


/ 

0.5 


u. 


a(z) 


dz 


JL • ^ 

/ a(z) dz 
0.5 


• • • » 


( 22 ) 


from the use of equation (12) and (16) the above relation can be 
reduced ±o, 


T) 



du 


1 

dz 


The yield is given by , 


• • • • 


(23) 


Y = Ti X S (24) 

Where Y is the yield of the desired product. The yield 
represents the actual amount of B formed and it may be possible 
to counteract low selectivity with high activity to maintain the 
yield. 

The differential equations (6) and (7) were solved , 
and it was found that at the dimensionless radius equal to 0.72 
the flux becomes zero. Keeping this' in view we define a general 
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function which satisfies equation (12). The function is given 
as, 

23-“x(|)3-“ _ 

9 ( 2 ) = X z “ (3-a) 

^3-a /4x3-cc 

"'‘3'' when z<0.72 

.... (25) 

z"^® (a+3) 

when 2 >_ 0.72 

where a ^ ^ 

When we change the values of alpha, we got the different 
type of distributions; 

a =0 uniform 

a = 1 linear 

a =2 parabolic etc. 

By using these expressions we solved the equations for 
various values of a and the results are given in chapter 4 and 
the method of solution is discussed in chapter 3. 

DIRAC DELTA TYPE DISTRIBUTION; 

After analysing the various distribution functions, 
now we will analyse the Dirac - delta type distribution. For 


37 

192 




1.0 
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the Amphora II, we define the function which must satisfy the 
constant of equation (12) : 


a(z) 


d (z-z^) 

24 2 

7 ^ 


. ... (26) 


Amphora shape has two surfaces, and the diffusion 
takes place from both the surfaces. Accordingly there will be 
two Dirac-delta functions, one for external surface and another 
for internal surface, stisfying the requirement of the charact- 
eristic equation (12): 

Substituting equation (26) in equation (6) and (7) 
we get the following relations. 



2^ ^ 
dz 

= “i; 

c 

(27) 

Or 

1.0 

/ dUj^ 

^c 


1.0 . 

/ . z dz 

^cext 



^‘"l.C 

= “l,c 

^ext 

...... (28) 


Where u, _ is the concentration of component A , where 
j. , c 

active sites are situated. And z^^^s the active catalyst location, 
and is given by , 


= / z“^ dz 

^©ext 


(29) 


and for the internal surface it becomes, 

^cint^ 

u.. . = f tT dz 

0.5 


( 30 ) 
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Similarly equation (7) becomes • 


1“ u. 


^ext ^^2 '^2,c 




( 31 ) 


where . is concentration of component B at the location 
where the active material is applied. 

Now we put this Dirac-delta function in equation (21) 
and we get the global selectivity for this function. 


S = 


/ u«(d(z-z ) / z^) dz 
( A ) 2C ^ ■ 

J z^ u, (d(z-z^) / z^)dz 


(32) 


S = 


(4) 


(33) 


And also the selectivity is defined by equation (18), 


s = - (4) 


(dU 2 / dz) 
(du,/dz) 


(18) 


equation (18) can be written as. 


s = - Ci) (557) 


(34) 


putting the values of derivetives the above equation can be 
reduced. 


__ i - u^o 

S = i-^) ^ 

1.0-ui,c 


( 35 ) 
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Equation (33) and (35) are the same, so equating these we have, 


( 4 ) = 1 . 0 - ( 4 ) 


-u 


1 * 0— u 


l,c 


l,c 


(36) 


Rearranging the above equation. 


*'2,c = “l.c < 




(37) 


In the above relation we have written u^ in terms of 

Zf c 

u, Thus a function G can be defined such as u^ =G(u, ) 

J. C Z j C X c 


Hence, 


•E-T^u, _ +1,0 

= “i,c > 


(38) 


equation (34) can be rewritten in the form of function G Thence 
the maximum selectivity is given by, 


G (u 




) - u. 


max 




(39) 


or equivalently. 


max 


1.0 - ( 4 -) 


G(u 


l,opt^ 


(40) 


u 


1, opt 


Where Uj^^p^ is the ^ value which maximizes the right hand 
sides of equation (39) and (40). From equations (36)- (40) we can 
determine the value of and Uj^ opt' graphically. The tengeni 

gives the value of Uj^ at the particular operating conditions. 
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Analitically the slope of the tangent can be given as 


- 


max* 


Since has been determined and now it is 

which is to be determine. For external and internal surfaces 
the optimal parameters of equations (28) and (29) are given, 


and 


^^opt^ext 


^^opt^int 


37 

192 


1 .0~u 


l,opt 


^1 '^IjOpt 


= d 


ext 
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“l,opt 


“ ^int 


(41) 


(42) 


Putting the values of the above optimized parameters into 
equations (28) and (29), we have, 


d 


ext ~ 


1.0 

/ z”^ dz 

^^opt^ext 


Or, 


^ ^opt^ext 


l.d+ d 


ext 


.... (43) 


and 


^int 


^^opt ^int 
/ z ”2 dz 

0.5 


Or, 


^^opt^int 


2-0 - ^int 


# # • • 


(44) 


Thus the optimal length of location of a thin layer 
inside the catalyst is given by equations (43) and (44) for « 
external and internal surfaces respectively. It is also clear 
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from the labove equations that as the thiele modulus increases, 
the optimal location shifts towards the surface, and at very 
high values of the thiele modulus, the whole layer shifts at 
the surface, thus the egg shell type distribution is the 
best at these conditions, 

STEP DISTRIBUTION ; 

Practically it may not be possible and also from 
catalyst dispersion and sintering viewpoints, to: locate all 
the active catalyst at one specific location, as was described 
in the case of Dirac - delta type distribution. We now examine 
a more realistic case in. which the catalyst is present only in 
a narrow region and centered, of course, about at the optimal 
distance as was previously determined for amphora. Now the 
distribution function is defined by the following equations. 



a(z) 

= ° ^ < ^opt 

z>Z2^ 

.... (46) 



= P ^opt 


.... (47) 

where 

^opt 

<1 

1 

In 

II 


(48) 

and 

^1 

= z +A 


.... (49) 

where 

A is the 

width of the step 

function. 



We define the p such as to satisfying the requirement 
of equation (12), Thus p comes out to be 

p = I (• ) ,... (50) 

^1 ” ^opt 
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SIMULTANEOUS REACTIONS; 

Consider the following case of parallel first 
order, isothermal reactions with negligible external mass 
transfer resistance. 



With the rate constants Kj, and K 2 * The steady state 
mass balance equations, in the nondimensional form, with the 
dimensionless variables as were defined in the case of 
consecutive reactions are given, 

dilj 

^ 3z dz^ ^ “1+ “2' 

y --csa) 

z dz 

And the non dimensional boundary conditions are given 


by the following relations, 

u^(l) = Uj^CO.S) = 1.0 (53) 

U2(l) = U2(0.5) =U2'^ (54) 


The global selectivity S is defined by equation (17) and (18). 
/^plication of the steady state mass balance equations into 
equation (17)and (18) gives the relation for the global selecti- 
vity for simultaneous reactions as follows: 
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1.0 « « 

/ 07 a(z) dz 

_ 0.5 

S = Q .... ( 55 ) 

/ * { •^ 0 ^ ^ 1+^2 ^ 2 ^ 

0.5 

The steady state mass balance equation with a(z) ( as 
defined by equation (25)) are solved by shooting method (described 
in chapter 3) and the selectivity for uniform and linear, parabolic 
distribution has been evolved by using equation (49). The results 
are given in chapter (4). 

The effectiveness factor for simultaneous reactions 
is given by their defining relationship which reduces to , 

^ ^ ( <^1 ^1+4 *^2 ^2 ) z^ a(z) dz .... (56) 

^1 0.5 ^ 

and the yield, 

Y = TixS (57) 

The results are discussed in chapter (4). 

DIRAC-DELTA TYPE DISTRIBUTION: 

The Dirac-delta type distribution is defined by 
equation (26) and the optimization technique is also the same as 
was described in the case of consecutive reactions. But the 
optimal parameters differ that of the consecutive case. The 
maximum selectivity is given by the following equation. 
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^ _ ^l,opt ^2 

ni3x 2 o 

'^l,opA ■‘■^2 ^2, opt 


(58) 


Where Uj^ is the value of which maximizes 

the selectivity function and hence the optimal Dirac distri- 
bution is given 
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T 



(59) 


and thus the maximum global selectivity can be obtained by a 
Dirac-delta type distribution and the optimal distance to apply 
the thin layer of the active sites is given, 


^^opt^ext 


1.0 


1+ d 


(60 ) 


ext 




1.0 


opt'int 2.0 -c 


(61 ) 


int 


where 

* 1.0 n 

^ext X 


39.0 

r^o 


(02 + 0|/^) 


... (62) 


* l^opt^int 

d._.u = / 2"2 dz 


int 


0.5 


17.0 (1« 0^1^ opt) 

192.0.^2 


(^1 + <^ 2 /"^) ^l,opt 


... (63) 


Again the egg shell type distribution is the best 
when thiele parameters are very high. 
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METHOD OF SOLUTION 

In formulation of the problem we got two system of 
differential equations, one for consecutive reaction » 
equation (6) and (7) along with the boundary conditions (8) 
and (9) and the other for simultaneous reactions, equations (51) 
and (52) along with the boundary conditions (53) and (54), These 
ere highly nonlinear equations, and the analytical solutions are 
not possible. We employed the numerical techniques to solve 
these problems, which provide fairly accurate results. 

Number of approximate methods are avilable for 
solution of differential equations. They can be classified 
under two categories, 

(1) Method of weighted residuals. 

(2) The variational method. 

In the method of weighted residuals one works directly 
with the differential equations and boundary conditions whereas in 
variational method one uses a functional related to differential 
equation and boundary conditions, i.e., the solution is expanded , 

I 

in terms of a trial function with undetermined constants or i 

functions. The method of weighted residuals provides a framework 
to compare and elucidate the features of individual method, which it 
encompasses e.g. collocation, Galarkin, integral etc. 
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Variational methods are not applicable to all 
problems and thus suffer a lack of generality. Sometimes 
they provide powerful results, such as upper and lower bounds 
on qualities of interest. e. g. , a complementry variational princi- 
ple method gives correct bounds to effectiveness factor. 

Further, methods of weighted residuals are easy to apply whereas 
variational methods require more mathematical manipulations. 

For the above reasons, the method of weighted residu- 
als is considered here, though all the approximate solutions 
satisfy only part of the conditions of the problem; for example 
the differential equation may be solved at few positions, rather 
than at each point. The method of weighted residuals encompases 
several methods as given above, i.e., collocation, Galarkin, 
finite difference, finite element, method of moments, the least 
squares etc. A brief description of these methods is given below. 

In the finite difference computations the difference 
expressions corresponds to a linear interpolation between grid 
points. When the solution at the grid points is known the 
solution at other positions is found by linear interpolation. 
Other more sophisticated schemes are possible. With spline 
functions, the region between two grid points is interpolated 
with second or higher degree polynomials. The derivatives at 
the grid points can be made continuous. Examples of these types 
of approximations are given by Hall and Watt (15) for integrating 
non linear ordinary: differential equation as. two point boundary 
value problem. The finite element method differs from finite 
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difference method and other methods, in that with spline 
functions a low order polynomial is used in each region from 
^i Least squares method minimizes the mean 

square residuals of the problem. Galarkin's method is simi- 
lar to collocation method but is very useful for instability 
or eigen value problems. Galarkin’s method requires the 
evaluation of integrals which could be difficult depending upon 
nature of integrand, though quadrature method can be applied 
to ease this problem. 

In the collocation technique, it is only necessary 
to evaluate the residuals at the collocation points, which 
could be done very easily even if a highly non linear term 
appears in the differential equation. This method differs in 
approach from other methods in that it uses a single higher 
order polynomial over the whole range of interest Xj^ ^ x ^ x^^, 
and this polynomial is arranged to satisfy the boundary conditions. 
This method can be used to initial value problems using Legendre 
polynomials as spline functions (13). The collocation technique 
was first advanced by Lanezos. Harvay and spiess (46), used 
polynomials which were orthogonal on the boundary. Villadsen 
and Stewart (47) made the major advancemerft; when they developed 
orthogonal collocation for boundary value problems. 

They choose the trial functions to be sets of ortho- 
gonal polynomials which satisfied the boundary conditions and 
the roots of the polynomials gave the collocation points. The 
orthogonal polynomials were considered because orthogonality 
takes care of uniqueness of roots and reality of the roots. The 
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choice of collocation points thus become fixed than arbitrary 
and the low order collocation results are more dependable. 

A major simplification can be done if solution is desired in 
terms of coefficients of trial function. The whole problem is 
then reduced to a set of matrix equations which are easily 
generated and solved on computer. Quadrature formulae are also 
provided which are especially important if the primary informa- 
tion desired from the solution is some integrated property. 

Some other specialized methods are also available 
for the problem of simultaneous diffusion and reaction in 
catalyst pores. A brief description of these are worth to 
mention. 


The well known quasilinearization technique requires 
the partial derivatives of right hand side terms of diffusion 
equation, i.e., reaction rate terms withrespect to each variable. 
For some cases this terms may be such that evaluation of partial 
derivatives become very difficult, for example in the case of 
Langmuir- Hinshelwood type of kinetics. Also quasilinearization 
technique is not suitable where sharp maximum occurs. In another 
method which is known as invariant embedding, missing condition, 
usually the slope, at the known point is calculated. The whole 
region is divided into elements and from the newly calculated 
first point, a next point is calculated and again the procedure 
of finding the first point is repeated. 

A second order approximation is used to findout the 
the next point. Agarwal and Saraf (1) have described this. They 
used the method in solving simultaneous two point boundary value 



problems. It is also shown that computer time doubles for 
each increase in:the number of division by one. 

CONSECUTIVE REACTIONS ; 

The non linear two point boundary value problems, as 
we obtained in chapter 2 are solved by finite difference method, 
(8,13,15,48). 


The application of finite difference method to boundary 
value problem is straight forward. Any derivative occuring in 
the differential equation are replaced by their finite difference 
approximation. The steady state mass balance equations, as were 
reported : 


1 


du 


dz 


(z^ ^)= 02 a(z) U 2 


( 6 ) 


can be reduced to 


d^u- 


dz' 


2 ill 

^ z dz 


02 a(z) U 2 =0 


(64) 


and similarly the mass balance equations for component B ean be 
reduced to , 

2 

d u_ du^ - « 

i z dz “^2 ^2 a(z) ^2 = O ... (65) 

dz 

It is clear from above equations that they are 
simultaneous. Hence from the equations (64), U 2 can be calculated 
at the specified values of 'z* and then at this specified value 
of z, U 2 can be calculated by putting the value of Uj^, in > 

equation (65). The finite difference approximation gives for i 



equation (64) , 
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hxzn 

- 0f a(z^)(u^)^ =0 (66) 

where the mesh is taken uniform and is defined by the following 


equation, 

Xn = a+nh, 0 _< n _< N (67) 

and h = ^ (68) 


where b and a are the two boundary points, which are given 
by equations (8) and (9). it is obvious that we got a band 
matrix of algabraic equations, which can be solved by a 
subroutine BAND and the value of (^ 2 )^^ ai'e calculated by 
the subroutine DIFFl. 

Equation (65) may be written in the form of approx- 
mate finite difference equations such as: 

02(ui)n a(z,i) = 0 (69) 

The above equation can be solved by using the subrou 
tine DIFF2 by calling BAND matrix subroutine. 
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After calculating the values of the and U 2 equa- 
tion (21) can be numerically integrated by putting the values 
of (u 2 )j^ and (u 2 )j^ at the difference points and thus the 

global selectivity can be calculated. 

For isolated solutions of a nonlinear problem the 
difference scheme has, for sufficiently small h, a unique 
solution converging to the exact solution if, 

(1) The linearized difference equations are stable and 
consistent for the linearized problem, 

(2) The linearized difference operator is lispschitz continuous 

(3) The nonlinear difference equations are consistent with 
the non linear differential equation. 

SOLUTION OF THE ALGABRIC EQUATIONS; 

If the differential equation and boundary conditions 
are linear , so are the equations of the discretization. The 
equations (66) and (69) gives the triple-diagonal set of 
equations, and the matrix of these equations display the 
typical structure obtained from finite difference methods. 

If a low order approximation is used, the matrix will have 
a narrow band width, but for good accuracy h must be small and 
the matrix will have large dimensions. The solution of the 
Band matrix can be obtained by the subroutine BAND. 

SIMULTANEOUS REACTIONS t 

The mass balance equations (51 and 52), obtained 
in the case of simultaneous reaction can be solved by 
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relaxation, invariant embodient and shooting method. Here 
we employed the shooting method to solve these boundary 
value problems. 


SHOOTING METHOD ; 

In the shooting method we choose values for all the 
dependent variables at one boundary. These values must be 
consistant with any lOther boundary conditions for that 
boundary, but otherwise are arranged to depend on arbitrary 
free parameters whose values we initially guess. We then 
integrate the odinary differential equations by initial value 
methods, arriving at the other boundary. And then iterate till 
the approximate solution is arrived. 

The steady state mass balance equations. 


du 




with the boundary conditions 


and 


Uj^(l) = Uj^(0.5) =1.0 

U 2 (l) = 02(0.5) =U 2 ° 


( 70 ) 

( 71 ) 


( 72 ) 

( 73 ) 


Can be reduced into four first order boundary 
value problems, such as. 




( 74 ) 


dU2 

dz 



( 75 ) 
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dUo 0 9 10 

+ I U 3 = (0J Ui+ i 0^ Uj) a(z) .... (76) 

dua 00 • 

^ + 7 ^ 4 =“ ^ 0 lUjLa(z)) .... (77) 

with boundary conditions 
Uj^(l) = u^(0.5) = 1.0 

U2(l) =U2(0.5) = U2° 

In shooting method we choose the initial guesses 
of the remaining variables u^and and define a vector 
which contain the boundary conditions at the inner surface. 
These boundary conditions are feed in the form of a subroutine 
ODIENT. 


- \ 

We define vector F^on the outer surface. The values 
are feeded in the form of a subroutine SCORE. As far as the 
Newton-Raphson is concerned, we are nearly in business. We 
want to find a vector value which zeros the vactor value of 

— V 

F . We do it by computing the solution of two linear equations, 


[e]dV^ - F^ .... (78) 

and then adding the correction back 

Hr^new ^ -^Id ^ (79) 


and the Jacobian is defined as 




..... (80) 


To solve the initial value problems fourth order 
Range Kutta method is employed . The description of which 
is given in standard texts (13,48). 
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RESULTS AND DISCUSSIONS 


In this section we will discuss the effect of various 
ph^siochemical parameters on selectivity, defined by equations (19) 
and (55), effectiveness factor, given by equations (22] and (56) and 
yield defined by equations (24) and (57), in the case of consecu- 
tive and simultaneous reactions. The range of different parameters 
studied are given in table (4) and are defined below, 

(1) THIELE MODULUS (0 ) 

It is defined as the ratio of reaction rate to diffusion 
rate. It represents the resistance to diffusion in the pellet. 
Higher the value of the thiele modulus, means, the reaction is 
severely diffusion limited. The results show that at higher 
value of this parameter the effectiveness factor and the selecti- 
vity are poor. 

(2) ACTIVITY DISTRIBUTION PARAMETER (a) 

This parameter gives the distribution of activity in a 
catalyst. By changing this parameter, we can get the different 
type of distribution, 

a = 0 uniform 

a = 1 linear 

a = 2 parabolic 

a > 2 other distribution 

The values of a . have been varied in the range of (D to 8. It 
is found that higher values of a gives higher value of selectivity 



and effectiveness factor 
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(3) RATIO OF KINETIC CONSTA^JTS rA^ 

was defined by equation (20). We studied the 
effect of ^ on selectivity, effectiveness factor and yield. 

It is observed that higher value of ^ gives poor selectivity, 
yield and effectiveness factor. 

(4) DIMENSIONLESS SURFACE CONCENTRATION OF COMPONENT B 

It is defined as the ratio of surface concentration 
of component B to the surface concentration of component A. we 
studied effectiveness factor, selectivity and yield of the 
consecutive reactions as a function of u^ . We shown that lesser 
the concentration of u^* higher is the selectivity and effectivene 
factor and yield. 

(5) RATIO OF DIFFUSION COEFFICIENTS ) 

It is defined as Del/De2, where Del and De2 are the 
effective dif fusivities of the components A and B respectively. 

The effect of has been studied in the range of 0.1 to 10.0 . 

It was found that at higher values of this parameter, selecti- 
vity and yield can be improved. 



Table - 4 


RANGE OF PARAAETERS STUDIED 


Parameter 

1 . 0 , Thiele modulus 

2. a. Activity distribution parameter 

3. ratio of kinetic constants 

4. dimensionless surface concentration 
of component B 

5. .ratio of diffusivities 


Range 


0 . 1-100 

0-8 

01 - 1.0 

0 - 1.0 


0.1 - 10.0 





4a 



IT) 

& 

iZ 


Globol selectivity vs. thiele modulus for series reaction on 
amphoro cotolyst. Discontinuous curve shows the s for sphere. 
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CONSECUTIVE REACTIONS t 


The effectiveness factor as a function of thiele 
parameter is given in figure (4), with a as a second parameter. 
The uniform activity distribution (a=0) gives the lower value 
of effectiveness factor, as the value of a is increased, the 
effectiveness factor increases for all values of thiele para- 
meter. When the thiele parameter is less than unity, effecti- 
veness factor is close to unity. At higher value of thiele 
parameter (0 >3), the effectiveness factor decreases linearly 
with thiele parameter. The results on amphora catalyst has been 
compared with a conventional spherical pellet, and it has been 
found that we can improve the effectiveness factor tremendously 
by using amphora shape, at higher value of thiele modulus. At 
0j^=lO, T] exhibited by an sphere is only 0.11, while its value 
is 0.4 for uniformly active amphora and 0.55 for non uniformly 
active amphora (a=8). 

Figure (5) gives the selectivity as a function of thiele 
modulus for different values of a. Keeping and ^ as consta- 

nts. It is clear from the figure that as thiele modulus increase 
i.e., the reaction is severely diffusion limited, the selectivity 
decreases. But it can be improved by increasing the value of 
a. It represents that when the diffusion is severe, the activity 
should be at the decreasing order towards the centre of the 
pellet. 

For the same value of U 2 °and*T^ but different value of 
A ( the ratio of kinetic constants and Kj_) the results are 
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Global selectivity vs. thiele modulus for series reaction on 
omphoro cotolyst. 
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plotted in figure (6). The value of selectivity decreases 
when the value of /\ is increased from 0.1 to 0.3. It is also 
observed that the parameter a becomes more effective at higher 
value of <l>. 

Fig. (7) shows the yield of component B in conse- 
cutive reactions as a function of thiele parameter for different 
values of a. ( The value of =1.0, X =1.0, and>v =0,3) and 
Fig. (8) for (u^ = 1.0, =1.0 and 0,5). These two figures 

show a marked reduction of yield for high value of thiele 
parameter. Increasing the value of a. increases the yield, as 
before. 

In Fig. (9) selectivity is plotted as function of (The 
dimensionless surface concentration of component B), for different 
values of a, at 0j^ =10,0, When the surface concentration is zero, 
selectivity is high and it continuously decreases, when the 
surface concentration increases. This effect is important as 
in actual reactors, as the reaction proceeds, the value of u^ 
increases, and the selectivity, in the reactor goes down. The 
effect of change of U 2 on yield is shown in Fig. (10). 

Fig. (11) and (12) shows the effect of changing the 
value of (The ratio of diffusion coefficients) on selectivity 
and yield respectively. With increase in the value of ^ the 
value of selectivity and yield increases, also increasing the 
value of a improves both selectivity and activity. 
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Fig. (13) gives the effect of varying the value of;(( the 
ratio of kinetic constants Kj^ and K 2 ) selectivity for given 
values 02=10.0,'^ =1.0 and u^ =0.5. The results shows that 
when the value of is increased, the selectivity decreases. 

Also at high values of a, the selectivity is high. 

SIMULTANEOUS REACTIONS 

The effect of changing thiele parameter and tt(the 
activity distribution parameter) on effectiveness factor 
is shown in Fig. 14 for u^ = 1.0 ,t^ =1 and > = 0.1. The figure 
shows that as the thiele parameter is increased the activity 
decreases and at higher values of thiele parameter (0 > 3) the 
effectiveness factor decreases proportionally with 0 . The 
value of effectiveness factor increases by increasing the 
value of a. 

Fig. 15 gives the selectivity as a function of thiele 
parameter for two value of a when U 2 = 0 = l.O, ?>= O.l • 

This figure shows that selectivity is nearly unity at low value 
of 0 but decreases by small amoant act higher values of 0, 

(Thele parameter) , which is obvious also since both the reactions 
are of first order Fig. 16 gives the yield for the above case. 

Fig. (17) gives the selectivity as a function of u^. 

At lower values of u^ * the selectivity is high and it continuous! 
decreases, when the value of u^ increases. 

DIRAC-DELTA TYPE DISTRIBUTION 

The global selectivity is much higher in Drac-delta, 
than the global selectivity obtained, with other catalyst 
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distributions such as unfirom and linear. This has been shown 
in Fig (18) that when (2 qp^= 0.7, =0.1, u® =1.0, -^=1.0) the 

selectivity vs thiele modulus curve gives about 98;^ selectivity. 
As is shown in Fig. (19) that distance z^, where the active 
components are applied, is a function of thiele parameter. If 
thiele parameter is higher than the specified value, the 
selectivity becomes zero. 

Fig. (19) shows the Zqp^ as a function of thiele 
parameters at ^=1.0,"^= 1.0, 0^=1. 0). It is shown that 

when the reaction is kinetic controlled, the whole active sites 
should be situated at the centre of the pellet, the 'egg yolk' 
type distribution. As the thiele modulus is high (0 > 6*0) the 
optimal strategy is to apply the sites over the surface. 

EFFECTS OF THE PHYSIOCHEMICAL PARAMETERS IN GENERAL 

EFFECT OF u® : 

The effect of on selectivity, effectibeness factor 

and yield on both the parallel and series reactions is much 

telling. In case of the series reactions it plays an important 

role. Initially if u^ is zero, in the case of series reactions , 

the selectivity is high and it can never be negative. Higher 

the value of u® , in any case, lower will be the selectivity, 

2 

effectiveness factor and yield of the reactions. So tt is 
suggested that initially there should be no surface concentration 
of the component B. The effect of component B is more telling 
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in case of series reactions, since component B further dissoci- 
ates into the other components. 

EFFECT OF ^ ; 

?vis the ratio of the squares of the thiele modu- 
lus of the reactions. When is constant, it reduces -to the 
ratio of the kinetic parameters K 2 and Kj^ . At lower value of 
this parameter we gbtained higher values of selectivity and 
yield. The effect of is more telling in the case of series 
reactions. When = 1, tends that and hence the 

selectivity becomes zero. At higher values the selectivity 
becomes negative. But it can be improved, by using the Dirac- 
delta type distribution. 

EFFECT OF ; 

is a dimensionless parameter defined as the ratio 
of the effective diffusivity of component A to the effective 
diffusivity of component B. At lower value of this parameter 
we got lower selectivity and the effectiveness factor. 

EFFECT OF a: 

is a parameter defining the activity distribution. 
When a=0 the distribution is uniform, a=l, the distribution 
is called linear, a=2, the distribution is said to be parabolic. 
The higher value of the a represents that the activity is 
highly decreasing towards the centre of the pellet. As a 











Fig. 19 -Variation of optimal location againt the thiele 
modulus Dirac -delta distribution. 
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increases, selectivity yield and effectiveness factor is being 
increased. It effects, highly the parallel reactions, as it 
is clear from the Figures (4)"!» (10), 

COAS>ARISION WITH THE SPHERE : 

We compared the effectiveness factor and global 
selectivity obtained on the new shape, with sphere. It was 
found that at lower value of the thiele modulus, at the same 
operating conditions, both the sphere and the amphora exhibit 
the same values of effectiveness factor. But as thiele modulus 
increases, the amphora shape exhibits tremendously improved 
values of the thiele modulus. Regarding the selectivity, 
amphora catalyst, produces always, better results, at any 
operating condition. The comparision is made in Fig. (4) and (5). 

COMP ARI SION BETWEEN CONSECUTIVE AND SIMULTANEOUS REACTIONS: 

The non uniformly active catalyst exhibit improved 
selectivity in both the cases . But the effectiveness factor, 
changes much faster with changing the activity distribution 
parameter in the case of simultaneous reaction than the case of 
consecutive reactions^ but the yield and the selectivity changes 
tremendously in the case of consecutive reactions. It is due 
to that the selectivity in case of the consecutive reactions, 
depends very much of the intermediate product, since this 
component, decomposes further, to produce the non desired produ~ 
ct. V/hile in the case of simultaneous reactions, it is not so. 
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CONCLUSIONS 


Amphora catalyst pellet gives better activity 
selectivity and yield for consecutive reactions in compa- 
rision to spherical catalyst pellet* (e.g. 0j^=lO.O, 11=0.4 
and S = 0.66 for amphora and t)= 0.11 and S= 0.49 for 
spherical pellet). It gives higher selectivity in comparision 
to the conventional spherical pellet due to small diffusion 
path and low residance time of intermediates in the pellet. 

Thus for consecutive reactions it plays an important role at 
high values of the thiele parameter. 

It has been found that the amphora type pellet gives 
better activity , selectivity and yield for consecutive reactions. 

As the value of the activity distribution parameter in increased 
for 0 to a high value, there is an increase in activity selec- 
tivity and yield . ( e.g. , at a =0, t) =0.4, S =0.72; at a = 1.0 
Ti = 0.42, S 0.74) 

Among the various types of distributions studied, the 
Dirac-deTta type distribution gives the maximum selectivity, for 
a specified value of thiele parameter, ( at 1^= ICkQ. , A= 0.3, 

U 2 = 1.0 and =3.0 Dirac - delta, B8y. selectivity in comparisioi 
to the 6iy of uniformly active amphora and 63yi in comparision to 
the linear distribution) . It gives higher values of selectivity 
since it exposes the active material to the optimum of the reac- 
tants, and thus the secondary reactions can be minimized. This 
optimal distance is a function of thiele parameter . As the 
value of the thiele parameter increases, the optimum distance 
shifts towards the surface and at very high values, the 



64 


the strategy is to apply all the active material at the surface, 
which is known as *egg shell* type distribution . 

RECOMMENDATION FOR THE FUTURE WORK; 

The present work can be extended to more realistic 
cases like external phase mass transfer resistance, non isother- 
mal diffusion and reaction problem, convex and concave type 
profiles, application of non uniformly active amphora to the 
reactors, deactivation and ofcourse for a general order kinetics. 
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C THIS PROGRAMME CALCULATES THE SELECTIVITY ^EFFECTIVEiiNESSiK 

C FACTOR AND YIELD FOR ISOTHERMAL m FIRST ORDER SERIES % 

C REACTIONS <A >B >C) WITH NEGLIGIBLE EXTERNAL t 

C PHASE RESISTANCE ON NONLINIFORMLY ACTIVE AMPHORA SHAPE* * 

C THE SECOND ORDER r TWO POINT BOUNDARY VALUE PROBLEM HAS >lc 

C BEEN SOLVED BY USING THE FINITE DIFFERENCE METHOD* Ht 


C ALPHA IS A COSTT* ALWAYS GREATER THAN ZEF<0 sK 
C ZI IS THE RATIO OF EFFECTIVE DIFFUSIVITIES<D©A/DeB> t 
C ALEM IS THE RATIO OF THE SQUARES OF THIELE PARAMETERS ^ 
C (PHI 2/PHI !fc 
C PHIl AND PHI2 ARE THIELE PARAMETERS OF COMPONENTS 'A' * 
C AND "B' RESPECTIVELY* * 


C >1? * if? !!c ^ sK )!c IK >K >K 3^: He He 3t: * >K >i« ^ 5K >!c ^ >lc >K # >|i >!« 5^ )i« )K :ic >!e $ $ * ;!c >K 5}: * 5K )K * 5K He * JloK >!c )!c 5K * 

DIMENSION Yl( 100)1- Y2< 100) fX( 100) i-AK 1001-4) ^A2( 100f4) 
COMMON/SAI /ALPHA 
COMMON/BBC/PHI 1 
C0MM0N/CCC/PHI2 
COMMON/ABC/ZI 
COMMON/GDR/Yl 
C0MMQN/FRG/Y2 
EXTERNAL FN 

REAL H ? ALPHA i- STRIPS » ZI ? ALEM 
OPEN ( UN I T:==20 ? FILE:=’ " F0R20 * DAT ' ) 

OPEN ( UNIT=23 •/ FILE-= ' F0R23 * DAT ' ) 

READ ( 20 y He ) XA ? YAl y XB y YBl . 

READ < 20 y He) N 
READ ( 20 y He ) Y A2 y YB2 
STRIPS==-N-H 
Hr==(XE-“XA) /STRIPS 
ALPHA™0*0 
1 ZI''^=1*0 

800 ALEH™0*1 

700 WRITE<23ylll) 

WRITE(23yl01) 

WRITE ( 23 y 90) YAl 
WRITE <23? 120) ALPHA 
WRITE(23yl03) 21 
WRITE ( 23 y 104) ALEM 
WRITE ( 23 y 200) 

PHIl ===0*1 

500 X(1)====XA-H"1 

PHI2==SQRT(PHIlHePHIlHcALEH) 

CALL DIFFl ( Y1 y N y A1 y X y YAl y YBl y H ) 

CALL DIFF2 < Y2 y N y A2 y X y YA2 y YB2 y H ) 

CALL SIMPS ( SUM lyHyNyXAyXByXyFNyYAlyYBlyYl) 

CALL SIMPS (SUM2yHyNyXAyXByXyFNyYA2yYB2yY2) 
CHe$HeHeHeHeHeHeH!^HeHeHeHeH:He5HHe'^HeHcHeHcHeHcHeHeHtHeHeHcHe^HeHeHcHeHeHeHeHeHeHeHeHeHeHeHeHeHeHtHeHeH:HtHeH:5^HeHeH:HcHe>S 
C SGLOB IS THE GLOBAL SELECTIVITY DEFINED BY EQN*18$21 He 

C ETA IS THE EFFECTIVENESS FACTOR DEFINED BY EQN*23 He 

C YLD IS THE YIELD OF COMOPONENT 'B' DEFINED BY EQN*24 He 

C)KHeHe:^H:He^*H:HeHeHeHeHe#HeHeHeHeHeHeHe>KHeH!HeHeH!HeHe*HeHeH:HeHcHeHcHe*H:H:HeHeHcHeHeHe*Hc*HeHeHeHeH:HeHeHeHeHeHeHe)K 
SGL0B=1 *0" < (ALEMHeSUM2)/(ZIHeSUMl) ) ^ 

ETA====(24*/7* )HeSUMl 



YLD"SGLOB«ETA 

TYPE * ? SGLOB s- ETA y YLD 

WRITE < 23 n06> PHI 1 » SGLOB y ETA y YLD 

IF-(PHI1 4LT4 I 4 O) GO TO 404 

IF<PHI1 4 LT 4 2*0) GO TO 402 

GO TO 403 

402 PHI 1==PHI 1 + 1*0 
GO TO 500 

403 PHI1=PHI1+240 

IF<PHI1 4 LE 4 20*0) GO TO 500 
GO TO 405 

404 PHI1=PHI1>!c1040 
GO TO 500 

805 ALEM=ALEM+042 

IF<ALEH ♦LE* 1*0) GO TO 700 
ZI“ZI+1*0 

IF<ZI ♦LE* 3*0) GO TO 800 

ALPHA--=ALPHA+1 4 

IF (ALPHA 4LE» 2*0) GO TO 1 

ALPHA™ALPHA+240 

IF (ALPHA *LE4 1040) GO TO 1 

806 FORMAT (10XyF843y5XyF946y5XyF946y5XyF94 6) 

101 F0RMAT(/20Xy23HTHE SELECTIVITY RESULTS//) 

103 F0RMAT<5Xy 'ZI”' yF946) 

104 FORMAT ( 5X y 'ALEM~' yF946/) 

200 FORMAT (SXy 60 ( IH- ) /12X y 'PHIl'y9Xy 'SGLOB" ylOXy 'ETA' yllXy 'YLD'/ 
1 lOXySOdH--) ) 

1 20 FORMAT ( 5X y ' ALPHA== ' y F9 * 6 ) 

90 F0RMAT(5Xy41HTHE INITIAL CONCENTRATION OF COMPONENT B™yF 543 ) 

1 1 1 FORMAT ( 20X y 25 ( IHsfc ) ) 

STOP 

END 

C THIS SUBROUTINE CALLS THE SUBROUTINE BAMD3 AND SOLVES !5c 

C THE SECOND ORDER BOUNDARY VALUE PROBLEMS* IT GIVES THE )lc 

C VALUE OF THE CONCENTRATION OF COMPONENT 'A'* He 

C :!K $ :{c ;|c $ 5K >!c if: !fc 5H >K )K f. f. ?!c >K ^ # >K * >k 3jc )ic $ He '4( $ )ic sK ^ jK >ic He 5ic f 5K )K >K He 5f£ $ 5l£ ^ * !|c $ # 

SUBROUT INE D I FF 1 ( Y 1 y N y A 1 y X y YA 1 y YBl y H ) 

COHHON/BBC/PHIl 

DIMENSION Al(100y4) pYKIOO) yXdOO) 

DO 20 I==:l.yN 

A1 ( I y 1 ) ==^:l * O-H/2 4 HCFI ( X ( X ) > 

A 1 ( I y 2 ) ^=="- ( 2 4 O-HHcHHsGl ( X ( I ) ) ) 

A1 ( I y 3 ) =^=1 4 O+H/2 4 HeFl ( X ( I ) ) 

Al(Iy4)-HHeHH(FRl(X(I) ) 

20 X(I + 1)==X(I)+H 

A1 ( 1 y 4 ) ™A1 d y 4 ) -"Al ( 1 y 1 ) #YA1 

Aldyl)“0 4 0 

A 1 ( N y 4 > A 1 ( N y 4 ) - A 1 ( N y 3 ) He Y B 1 

Al(Ny3)-«0 40 

CALL BAND3(NyAlyYl) 

RETURN 

END 



c this subroutine calls subroutine bands and takes the * 

C VALES OF THE CONCENTRATION OF COMPONENT A FROM THOSE 5K 

C CALCULATED BY SUBROUTINE YAPl AND THEN SOLVES THE 

C BOUNDARY VALUE PROBLEM FOR THE CONCENTRAYION OF THE COM* 

C -PONENT 'B', * 

C*************************************************************** 
SUBROUTINE DIFF2 ( Y2 ? N ? A2 s- X ? YA2 y YB2 ? H ) 

COMMON/BEC/PHI:L 

C0MM0N/CCC/PHI2 

COMMON/GDR/Yl 

DIMENSION A2(100s-4> i-Y2( 100) yX< 100) ? Y1 ( 100) 

DO 10 I==lyN 

A2 ( I y 1 ) “1 « O-H/2 ♦ JlcFl ( X < I ) ) 

A2(I i-2)--~<2.0-H*H*G2<X(I)) ) 

A2(I?3)=l«0+H/2»0*Fl<xa)) 

A2 a ? 4 ) ==H*H*Y 1 ( I ) 5KFR2 < X ( I ) ) 

X(IH)™X<I)-fH 
10 CONTINUE 

A2a ?4)===A2<ly4)“ A2<l»l)*YA2 
A2<lyl)=0*0 

A2(Ny4)=A2CNy4)--A2<N»3)*YB2 

A2<N?3)~0*0 

CALL BAND3(NyA2yY2) 

RETURN 

END 

C**************************************************************** 
C**************************************************************** 
SUBROUTINE BANDSCNyAyX) 

DIMENSION A <100^4) 

DIMENSION X(IOO) 

H"N”"1 

DO 10 I™2yM 

A<Iy2)™A(Iy2)*A<I-ly2)-“A<Iyl)*ACI“-ly3) 

A<Iy3)=A<Iy3)*A<I~ly2) 

10 A(Iy4)=«A(Iy4)*A<I“lp2)“A<Iy l)*ACI“ly4) 

X(N) = (ACN?4))!cA(N“"l y2) -A(N?l )*A(N"-ly4))/<A(Ny2)* 

1 A ( N-l p 2 ) -A ( N y 1 ) *A ( N“1 y 3 ) ) 

DO 20 I==lyM 
J™N""I 

20 X( J) -"<AUy4)- A( Jy3)*X<J+l))/AC Jy2) 

RETURN 

END 

************************************************************** 

C THESE ARE THE GENERALISED FUNCTIONS WHICH ARE BEING * 

C USED IN THE SUBROUTINES YAPl AND YAP2* * 

C********************************************2J55lc**************** 
C*****************************************************^****5lc*** 

FUNCTION FKX) 

Fl™2.0/X 

RETURN 

END 


C***********************************^************************** 



FUNCTION GKX) 

COHMON/BBC/PHIl 
COMMON/SAI /ALPHA 
PARA 1 ==2 0"" ALPHA ) 

PARA2“ ( A . 0/3 » 0 ) JKsk < 3 ♦ 0“ALPH A ) 

P ARA3== C A 0 0/3 ♦ 0 ) >K)Sc < 3 ♦ O-f ALPHA ) 

IFCX ♦LE* 0*725) GO TO 20 
GO TO 30 

20 G1““PHI i>KPHI 1* < 19.0/192* ):«: (3* 0-ALPHA > *PARAl5lsPARA2>k 

1 ( X)K^c (-ALPHA ) )/<PARAl -PAR A2) 

GO TO 40 

30 G 1“=“PHI iJfcPH 1 1>{C < 37 * 0/ 1 92 * ) >lc ( 3 * 0+ ALPHA ) !KPARA3)K ( X)!«)KALPHA ) 

1 /(PARA3-1.0) 

40 RETURN 

END 

C3RNc5ic*>l()K5KH«)lc He >!($*)!« He )i£)K>K#>l«!K5lt!<c5lc!!t^c)K)ic:)c#$:!c>!c)jc>lC!(C)K)!«)i{5lC)|cH:)tc:jc)i{![t)iJ 

FUNCTION G2(X) 

C0MM0N/CCC/PHI2 

COMMON/SAI/ALPHA 

PARA 1-2 * OHcHe (3 * 0- ALPHA ) 

PARA2- ( 4 . 0/3 . 0 ) H!He < 3 . 0-ALPHA ) 

P ARA3- < 4 * 0/3 . 0 ) *He ( 3 * 0-fALPHA > 

IF(X *LE* 0*725) GO TO 20 
GO TO 30 

20 G2™-PHI2HePH I2>K ( 19 * 0/ 192 * ) He ( 3 * 0-ALPHA ) HePARAlHePARA2He 

1 (XHeHe(-ALPHA) )/(PARAl-PARA2) 

GO TO 40 

30 G2--PHI 2HePH I2He ( 37 . 0/ 1 92 ♦ 0 ) He < 3 * O FALPH A > HePARA3He < XHeHc 

1 ALPHA)/CPARA3-1 *0) 

40 RETURN 

END 

CHeHeHeHeHeHeHeHeHe>KHeHeHeHeHeHeHeHeHeHe^^HeHeHeHeHeHeHeHeHcHeHe>KHeHeHe>K^$HeHeHtH!HeHeHcHeHeHeHeHeHeH:H:HeHeHe*He 

CHe H< He He He He He He He He He He H; He He He He )K He He He He He He He He He H; :’te He He He He He He He H: He He 5K He He He >K He He He sK He He >F. He He He He He He He H: He 
FUNCTION FRKX) 

FR1='=0*0 

RETURN 

END 

CHe H: He He He 3ie He He He He He He He He He He He 5^^He He He He He He He He He He He He He He He He He He He He He He He He He*H:He He He He He He He He He He 

CHe:ieHe;HeHeHeHeHeHeHeHeHeHeHeHeH5KHeHeHe:^.HeHeHe5KHeHeHeHeHeHeHeHeHeHeHesKHeHe)KHeHeHeHeHH-HeHeHeHe>KHeHeHeHeHeHeHeHeHe 
FUNCTION FR2<X) 

COMHON/SAI/ALPHA 
COHMON/BBC/PHIl 
COMMON/ ABC/Z I 
P AR A 1 -2 * one He ( 3 * 0- ALPHA ) 

P AR A2== ( 4 * 0/3 * 0 ) HeHe ( 3 * 0-ALPHA ) 

PARA3- ( 4 * 0/3 * 0 ) *He ( 3 * OF ALPHA ) 

IF(X *LE* 0*725) GO TO 20 
GO TO 30 

20 FR2--PH 1 1 HePH I IHeZ I He < 1 9 * 0/ 1 92 * 0 ) He < 3 * 0-ALPH A ) HeP AR A1 

1 HePARA2He(XHeHe(-ALPHA) )/(PARAl-PARA2) 

GO TO 40 

30 FR2--==-PHI IHePHI lHeZIHe(37 * 0/192 * 0 ) He (3 * OFALPHA )HePARA3He 

1 ( XHeHe ALPHA ) / < PARA3-1 * 0 ) 

40 RETURN 



END 

C THIS SUBROUTINE INTEGRATES THE EQUATIONS t 

C BY USING THE SIMPSONS' COMPOSITE RULE. :ic 

C>ic:^#)!c5!«5K>ic*)!{St5if>i«*$:lC5S<)loK5i«*>K«!K*^*3|cH5!f{5K>H5K>K>i?H«>i!^)}c)!c:|{5K5it:it^*3|C5k#H!5K*>ls!}c>ic>i«>ic*Mc 
SUBROUT I NE S I MPS < SUM ? H ? N y XA y XB i- X » FN i- Y A r YB i- Y ) 

DIMENSION XCiOO) vYdOO) 

ri™fn(xa:>)|cya 

R2==FN<XB)*YB 

H0V2“H/2.0 

HALF”FN<XA-fH0M2) 

S=0.0 
DO 10 

S“S4-Y(I)5{cFN(X<I)) 

HALF”H ALF+Y < I ) *FN < X < I) +1-1002 ) 

X(I + 1):=X(I>+H 
10 CONTINUE 

SUM® ( H/6 . 0 ) * < R 1 •fR2+4 . 0:f:HALF+2 . 0*S > 

RETURN 

END 

C ACTIVITY DISTRIBUTION FUNCTION A(X> t 

FUNCTION A(X) 

COMMON/SAI/ALPHA 
PARA 1 ®2 . Om ( 3 . 0- ALPHA ) 

PARA2® ( 4 . 0/3 . 0 ) >ic5!c <3 . 0-ALPHA ) 

P ARA3® ( 4 . 0/3 . 0 )!!««< 3 . 0+ALPHA ) 

IF(X .LE. 0.725) GO TO 20 
GO TO 30 

20 A® ( 1 9 . 0/ 192 . ) 5K ( 3 . 0-ALPH A ) JRP ARAl*PARA2>i{ ( XJifJlc ( -ALPHA > ) 

1 /(PARA1™PARA2> 

GO TO 40 

30 A® ( 37 . 0/192 . ) )ic < 3 . 0+ALPHA ) JKP ARA3)ic ( X>i£:lcALPHA ) / ( P ARA3-1 . 0 ) 

40 RETURN 

END 

FUNCTION FN<X) 

FN®A<X):^cX3i«X 

RETURN 

END 



:alculates the selectivity s-EFFECTIMEHESS factor * 

C AND YIELD FOR ISOTHERMAL? FIRST ORDER REACTIONS A 

C A >C? WITH NEGLIGIBLE EXTERNAL PHASE MASS TRANSFER 

C RESISTENCE* V<1) AND V<2) ARE THE MISSING BOUNDARY CONDITIONS 

C AND ALSO THE CONCENTRATION GRADIENTS ON THE SURFACE^THE REMAINING 

C NOMENCLATURE IS THE SAME AS WE UTILISED IN THE LAST PROGRAMME* 

DIMENSION V<2)?DV<2) 

COMMON/SAI/ALPHA 
COMMON/BBC/PHI 1 ? PHI2 ? ZI 
OPEN <UNIT=:23? FILE- 'YY2* OUT') 

YA2==1*0 
ALPHA=0*0 
ZI-1.0 
ALEM=0 * 1 
WRITE(23? Ill) 

WRITE(23?101) 

WRITE(23?90)YA2 
WRITE(23?120) ALPHA 
WRITE <23? 103) ZI 
WRITE <23? 200) 

PHIi=»0*l 

505 PHI2-SQRT <PHIl!icPHI l:icALEM) 

CALL SIMUL<y?DV?EPS?Xl?X2?Hl?HMIN?N) 

SGL0B«“V<2)/<V<1))KZI) 

ET A- < 24 ♦ 0/7 ♦ 0 ) *V < 1 ) / < PHI 1 HcPHI 1 •1PHI2*PH 12 ) 

YLD==SGLOB>i{ETA 

WRITE < * ? )!? ) PHI 1 ? SGLOB ? ETA 

WRITE <23? 106) PHIl ? SGLOB? ETA? YLD 

IF<PHI1*LT*1*0) GO TO 404 

IF<PHI1*LT*2*0) GO TO 402 

GO TO 403 

402 PH II “PH II +1*0 
GO TO 505 

403 PHIl-=PHIl-f2*0 
IF<PHI1*LE*20*0) GO TO 505 
GO TO 405 

404 PHIl“PHIl)lcl0*0 
GO TO 505 

405 STOP 

106 F0RMAT<10X?F8.3?5X?F9*6?5X?F9*6?5X?F9*6) 

101 FORMAT </20X?39HTHE SELECTIVITY ON AMPHORA J PARALLEL RXN//) 

103 F0RMAT<5X? 'ZI=' ?F9*6) 

104 F0RMAT<5X?'ALEM-'?F9*6) 


FORMAT < 5X ? 60 < 1H-) /12X ? ' PHI 1 


'SGLOB' ?10X? 'ETA' ?11X? 'YLD'/ 


1 10X?50<1H")) 

120 F0RMAT<5X? 'ALPHA-' ?F9*6) 

90 F0RMAT<5X?41HTHE INITIAL CONCENTRATION OF COMPONENT B-?F5.3) 

111 F0RMAT<20X?25<1H>K)) 

END 

C THIS SUBROUTINE CALCULATES THE MISSING BOUNDARY CONDITIONS He 

C V<1) AND V<2) AND RETURNS THEM TO THE MAIN PROGRAMME* He 

C N2“NUMBER OF MISSING BOUNDARY CONDITIONS He 


^ ^ ^ 



C DELTA IS THE PERMISSIBLE ERROR IN THE CALCULATED UALUES He 
C*HcHc*HcHcHc>fc*HcHc)k5KHc«>KHeHcHc*)ieH{*)kHc*He>!cH«)leH<He****TO>!c)K)K3icHc>|{3k>K)|{HcHcHc>KHc5iC5!{)lc)KJK)ie)ieH{)lc>it)5cHc 
SUBROUTINE SIMUL(Vi-DV?EPSfXl ?X2!-H1 H-IMIN^N) 

PARAMETER ( NyAR=4 . N2=2 » DELTA® 1 * OE-3 y EPS® 1 » 0E~6 y DX® 1 ♦ OE-4 ) 
DIMENSION V(2)yriELV(2>yF(2)yDU<2)yY<4) 

OPEN < UNIT®23 y FI LE® ' F0R23 ♦ DAT ' ) 

V(l)®0*4 

V(2)®0.4 

DELV(1)®DELTA5|cV(1) 

DELV<2>®DELV(1) 

Xl®0«5 

X2®1«0 

Hl®0*025 

HMIN=0*0 

N®4 

2 CALL SH00T<NVARyVyDELVyN2yXlyX2yEPSyHlyHMINyFyDV) 

IF(<ABS<DV<1.) ) ♦GT^ABSa-PSHeUa) ) ) ♦OR* 

* (<DU(2))*GT*AES(EPS>KV(2>))) GOTO 2 

RETURN 
END 

C THIS SUBROUTINE CONTAINS THE KNOWN BC'S AT THE INNER SURFACE 

SUBROUTINE LOAD<X:L yU y Y) 

DIMENSION U<2)yY<4) 

Y(1)®1*0 
Y(2)®:U0 
Y<3)®V< 1) 

Y<4)®V(2) 

RETURN 

END 

C>KHcHcH<)KH:)i«)ic5|eHe>l£*>!e)!c)lcHe>KH{)l<>K>!c)lc>le>|c)i{)lC5icHeHc>|C)l(Hc>!«)lc)K)KJH5|tHc>K!lc!^)|£*Ht)!:)|<>K>K5|eHc>KH()|cH:#H:>|c>K**>lcHtH!He)|eH<* 

C THIS SUBROUTINE CONTAINS THE BC'S AT THE OUTER SURFACE 
C>leHcHe>KHcHeHcHeHc)!c)!cH:He^HeHeHe>lcHeHcHeHe>ic!^Hc)l<Hc>!!)K)K)K)KHeHcH!HeHeHcHtHC5i(5ltHeHe5}cHt'^>lcHc)feHtH!Hc)K)le)KHeH!>le>KH!5l!He*Hc)leH::^ 
SUBROUTINE SC0RE<X2y YyF) 

DIMENSION Y(4)yF(2) 

F(1)®Y<3) 

F(2)®Y(4) 

RETURN 

END 

SUBROUTINE DERIVSCXy Y yDYDX) 

COMMON/BBC/PHI lyPHI2yZI 
DIMENSION Y<4>yDYDX(4) 

DYDX(1)®Y(3> 

DYDX(2)®Y(4) 

D YDX ( 3 ) ® < PH 1 2>KPH 1 25ieY ( 2 ) /Z I +PH 1 1 )1CPHI 1 « Y < 1 > ) H: ACT Y ( X ) -2 ♦ OH: Y ( 3 > /X 
DYDX < 4 ) ®“2 ♦ 0H:Y < 4 > /X~Y ( 1 ) HiZIsfePH I2HcPHI2H:ACTY ( X ) 

RETURN 

END 

SUBROUTINE SHOOT < NU AR y V y DELV y N2 y X 1 y X2 y EPS y H 1 y HM IN y F y D V ) 

PARAMETER (NP®20) 

DIMENSION V<N2) yDELV(N2) yF<N2> yDV(N2) yY<NP) yDFDV<NPyNP) y INDX<NP) 



EXTERNAL DERIVSfRKGC 
CALL L0AD(Xl5-g,Y) 

LALL SCORE (X2yY?F) 

»0 12 IV“lyN2 
SAV:=V(IV) 

V<IV)“V(IV)-MiELV<XU) 

CALL LOADCXl ?Mi»Y) 

CALL ODEINFCY yNVARyXl i'X2!'EPS!fHl yHMINsNOKfNBADfDERlVSyRKQC) 
CALL SCORE <X2yYyDV) 
no 11 1=1 yN2 

DFnv<iyiv>=(ny(i)-"F<i) )/nELV(iv) 

11 CONTINUE 
v(iv)=sAy 

12 CONTINUE 

DO 13 IV=lyN2 
DV<IV)=--F(IM) 

13 CONTINUE 

CALL LUDCMP ( DFDV y N2 y NP y INDX y DET ) 

CALL LUBKSB ( DFDU y N2 y NP y INDX y W ) 

DO 14 IV=lyN2 

V(IV)=V(IV)+DV<IV) 

14 CONTINUE 
RETURN 
END 

SUBROUT I NE RK4 < Y y DYDX y N y X y H y YOUT y DERI VS ) 

PARAMETER (NMAX=10) 

D IMENS ION Y < N ) y DYDX ( N > y YOUT < N ) y YT < NMAX ) y DYT < NMAX ) y D YM < NMAX > 

HH=H)1«0,5 

H6=H/6* 

XH=X+HH 
DO 11 1=1 yN 

YT < I ) =Y ( I ) THHjUDYDX < I ) 

11 CONTINUE 

CALL DERIVS<XHyYTyDYT) 

DO 12 1=1 yN 

YT<I)=Y(I)+HH>}cnYT(I) 

12 CONTINUE 

CALL DERIVS<XHyYTyDYM> 

DO 13 I=lyN 

YT(I)=Y(I)TH)KnYM<I) 

nYM<I)=DYT(I)+nYM(I) 

13 CONTINUE 

CALL nERIVS<X+HyYTyDYT) 

DO 14 I=lyN 

YOUT ( I ) =Y < I ) +H63IC ( DYDX < I ) +DYT < I ) +2 , )KDYM < I > ) 

14 CONTINUE 
RETURN 
END 

SUBROUTINE LUDCMP< AyNyNPyINDXyD) 

PARAMETER <NMAX=100yTINY=l *0E~20) 



DIMENSION A<NPyNP) » INDXCN ) f VV(NMAX) 

D==l* 

DO 12 I~1»N 
AAMAX~0. 

DO 11 J=1?N 

IF ( ABS<A( I p J) ) *01* AAMAX) AAMAX=ABS<A( I p J) ) 

11 CONTINUE 

IF (AAMAX*EQ»0*) PAUSE 'Singular matrij-u' 
VUCI)=1*/AAMAX 

12 CONTINUE 

DO 19 J=1pN 

IF <J*GT*1) THEN 
DO 14 I = 1pJ--1 
SUM==A<Ip J) 

IF (I*GT*1>THEN 
DO 13 K-1pI~1 

SUM=SUM“A ( I p K ) *A ( K p J ) 

13 CONTINUE 
A(IpJ)=SUM 

END IF 

14 CONTINUE 
ENDIF 
AAMAX~0* 

DO 16 I"JpN 
SUM«A(IpJ) 

IF (J^GTpDTHEN 
DO 15 K==1pJ"1 

SUM==SUM~A < I p K > 5fcA < K p J ) 

15 CONTINUE 
A(Ip J)™SUM 

ENDIF 

DUM“W(I))KABS<SUM) 

IF (DUM»GE,AAMAX) THEN 
IMAX==I 
AAMAX==DUM 
ENDIF 

16 CONTINUE 

IF < J*NE*IMAX)THEN 
DO 17 K™1pN 
DUM=A<IMAXpK> 

A<IMAXpK)"A<JpK> 

A( JpK)=DUM 

17 CONTINUE 
D=“D 

W<IMAX)=W( J> 

ENDIF 

INDXC J)--IMAX 
IF<J*NEpN)THEN 

IF(A<JpJ> ♦EQ*0*)A<J?J)=*TINY 
DUM-1*/A( Jp J) 

DO 18 I=:J-flpN 

A(Ip J)™A<IpJ)5KDUM 

18 CONTINUE 
ENDIF 

19 CONTINUE 



IF(A(N»N) ♦EQ*Ot ) A(Ni-N)-TINY 

RETURN 

END 

C LOWER DECOMPOSITION OF THE MATRICES 

C5|c!*5K*)l:!lC)ic*^*)|t)KHc)fc>K!t:*)!c)!(*****>|{*#J|{5K*JK>K5}:>K!jcJt!)K!K*^:5tC3l?Jt:!|?3!C!^5K*>i{«!C)^5!c)lt>f(5{c*#)K)f:Hc*)lc5K# 
SUBROUTINE LUBKSB ( A i- N v NP y INDX sr B ) 

DIMENSION A(NPyNP) »INDX(N) yB<N) 

11=0 

DO 12 I=lyN 
LL=INDX<I) 

SUM=B(LL) 

B(LL>=B<I> 

IF <II,NE.O>THEN 
DO 11 J=IIyI“l 

SUM=SUM“-A(IyJ)#B(J) 

11 CONTINUE 

ELSE IF (SUM*NE*0*) THEN 
II = I 
ENDIF 
B(I)=SUH 

12 CONTINUE 

DO 14 I=Nyly~l 
SUM=B<I) 

IF(I*LT»N)THEN 
DO 13 J=I+lyN 

SUM=SUM-"A(If J>*B(J> 

13 CONTINUE 
ENDIF 

B<I)=SUM/A<IyI) 

14 CONTINUE 
RETURN 
END 

SUBROUTINE ODEINT < YSTART y NUAR y XI y X2 y EPS y HI y HMIN y NOK y NB AD y DERI VS y RK 
>KQC) 

PARAMETER ( MAXSTP=10000 y NMAX=10 y TU0=2 » 0 y ZERQ=0 ♦ 0 y TINY=1 * E~30 ) 
COMMON /PATH/ KMAX y KOUNT yDXSAV yXP (200) y YP < lOy 200 > 

DIMENSION YSTART < NVAR > y YSCAL < NMAX ) y Y < NMAX ) y DYDX < NMAX > 

KMAX=20 

DXSAV=l*0E“-3 

X=X1 

H=SI0N(HlyX2“Xl) 

N0K=0 

NBAD=0 

K0UNT=0 

DO 11 I=lyNVAR 
Ya>=YSTART(I) 

11 CONTINUE 

XSAV=X~DXSAV*TWO 
DO 16 NSTP=lyMAXSTP 
CALL DERIVSCXyYyDYDX) 

DO 12 1=1 y NVAR 

YSCAL ( 1 ) =ABS ( Y ( I ) > TABS ( HjKDYDX < I ) ) +T INY 



12 CONTINUE 
1F(KMAX«8T*0)THEN 

I F < ABS < X-XSA V ) ♦ GT * ABS ( DXSAM ) > THEN 
IF ( KOUNT ♦ LT . KMAX-1 ) THEN 
KOUNT=KOUNT+1 
XP< KOUNT )”X 
DO 13 I-'IjNVAF? 

YP(IyKOUNT)~Y<I) 

13 CONTINUE 
XSAM“X 

ENDIF 

ENDIF 

ENDIF 

IF( <X+H“"X2))K(X+H™X1) ♦GT*ZERO> H=X2“X 

CALL RKQC < Y i- DYDX y NVAR » X H-l » EPS s- YSCAL sr HD ID HNEXT . DERIMS > 
IF<HDID*EQ*H>THEN 
NOK=NOK41 
ELSE 

NBAD==--NBAD+1 

ENDIF 

IF ( < X~X2 ) < X2~X1 ) ♦ GE * ZEr.RO ) THEN 
DO 14 I^l^NVAR 
YSTART<I>«=Y<I) 

14 CONTINUE 
IF(KMAX«NE*0)THEN 

KOUNT~KOUNT+l 
XP< KOUNT )”X 
DO 15 I=1»NMAR 
YP(I»KOUNT)"Y<I) 

15 CONTINUE 
ENDIF 
RETURN 

ENDIF 

IF< ABS<HNEXT) ♦LT«HMIN) PAUSE 'Stepisise smaller than miri;lmum+'' 
H-HNEXT 

16 CONTINUE 

PAUSE ''Too many steps*' 

RETURN 

END 

SUBROUTINE RKQC < Y s DYDX vUfKv HTR Y » EPS » YSC AL ^ HDI D » HNEXT » DERI VS ) 
PARAMETER (NMAX=10 i.FCOR= *0666666667 8- 
:<{ ONE-1* »SAFETY=0*9rERRC0N~6*E“4) 

EXTERNAL DERIVS 

D I MENSI ON Y < N ) r DYDX ( N ) i- YSCAL < N ) y YTEMP < NMAX ) i- YS AM ( NM AX ) » DYSAM < NMAX ) 
PGR0U™“0*20 
PSHRNK="0*25 
XSAM“X 
DO 11 I==l!-N 
YSAMCI)“Y<I) 

DYSAM(I>~DYDX<I) 

11 CONTINUE 
H“HTRY 
HH~0*5)KH 


1 



CALL RK4 ( YSAV y DYSAM i> N f XSAV » HH r YTEMP » DERI VS ) 

X==XSAV+HH 

CALL DERIVSCXi-YTEMPfDYDX) 

CALL RK4 < YTEMP y D YDX y N y X y HH y Y y DERI VS ) 

X“XSAV+H 

IF(X*EQ*XSAV)PAUSE ''Stepsis:© not si^niPicsnt in RKQC*'^ 

CALL RK4 ( YSAV y DYSAV y N y XSAV y H y YTEMP y DERI VS ) 

ERRMAX™0 ♦ 

DO 12 I=lyN 

YTEMP ( I ) ssY < I ) “YTEMP ( T ) 

ERRMAX“MAX < ERRMAX y ABS ( YTEMP ( I ) /YSCAL ( I ) > > 

12 CONTINUE 
ERRMAX=ERRMAX/EPS 

IF (ERRMAX ♦ST ♦ONE) THEN 

H“8 AFETYHcH* ( ERRMAX)!t'4!p3HRNK ) 

GOTO 1 
ELSE 
HDID“H 

I F ( ERRMAX ♦ GT ♦ ERRCON ) THEN 

HNEXT=SAFETYJHH)|c ( ERRMAX>!c>KPGROU ) 

ELSE 

HNEXT==4**H 
ENDIF 
END IF 

DO 13 I“lyN 

Y ( 1 ) =Y ( I > +YTEMP < I ) tfCQR 

13 CONTINUE 
RETURN 
END 

FUNCTION ACTY(X) 

COMMON/SAI/ALPHA 
PARAl =2 ♦ Om < 3 ♦ 0" ALPHA ) 

PAR A2= ( 4 ♦ 0/3 ♦ 0 ) 5K>1! < 3 ♦ 0~ ALPHA ) 

PAR A3== < 4 ♦ 0/3 ♦ 0+ ALPHA ) 

IF(X*LE*0»725) GO TO 20 
GO TO 30 

ACTY“ ( 19 ♦ 0/192 ♦ 0 ) >K ( 3 ♦ 0-ALPHA > '4cP ARA1JKPARA2* < (--ALPHA ) ) 

/(PARA1-PARA2) 

GO TO 40 

ACTY== ( 37 ♦ 0/192 ♦ 0 ) # ( 3 ♦ 0+ALPHA ) *P ARA3)K ( X#3KALPHA ) / ( PARA3--1 ♦ 0 ) 
RETURN 
END 


20 


30 

40 
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